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ABSTRACT 

Context. For accurately measuring intensities and determining magnetic field strengths of small-scale solar (magnetic) structure, 
knowledge of and compensation for the point spread function is crucial. For images recorded with the Swedish 1-meter Solar 
Telescope (SST), restoration with Multi-Frame Blind Deconvolution (MFBD) and Joint Phase Diverse Speckle (JPDS) methods 
lead to remarkable improvements in image quality but granulation contrasts that are too low, indicating additional stray light. 
Aims. We propose a method to compensate for stray light from high-order atmospheric aberrations not included in MFBD and JPDS 
processing. 

Methods. To compensate for uncorrected aberrations, a reformulation of the image restoration process is proposed that allows the 
average effect of hundreds of high-order modes to be compensated for by relying on Kolmogorov statistics for these modes. The 
applicability of the method requires simultaneous measurements of Fried's parameter r^. The method is tested with simulations as 
well as real data and extended to include compensation for conventional stray light. 

Results. We find that only part of the reduction of granulation contrast in SST images is due to uncompensated high-order aberrations. 
The remainder is still unaccounted for and attributed to stray light from the atmosphere, the telescope with its re-imaging system and 
to various high-altitude seeing effects. 

Conclusions. We conclude that statistical compensation of high-order modes is a viable method to reduce the loss of contrast oc- 
curring when a limited number of aberrations is explicitly compensated for with MFBD and JPDS processing. We show that good 
such compensation is possible with only 10 recorded frames. The main limitation of the method is that already MFBD and JPDS 
processing introduces high-order compensation that, if not taken into account, can lead to over-compensation. 
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and including stray light and the effects o f the large central 
Since the emergence of 3D simulations of solar convection, there obscur ation and spider, was determined by IWedemever-Bohml 
has remained a disturbing discrepancy between the measured d2008h . The obtained PSF leads to reconstructed granulation 
contrast of granulation and that obtained from simulations. It contrasts that are remarkably close to those of 3D simulations 
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, has been suspected that the major reason for this discrepancy dWedemever-Bohm & Rouppe van der Voortl I2009I) . A simpler 

O . is atmospheric and telescope stray light. However, difficulties fit of the PSF in terms of 4 Gaussians also leads to good 

' of accurately characterizing such stray light, in particular as re- agreements between granulati on contrast measured with BFI on 

^. \ gards the far wings of the corresponding point spread functions, Hinode and MHD simulations dMathew et alj2009l) . In addition, 

have prevented firm conclusions. Questions about the accuracy agreement between simulated images degraded to the resolu- 

of the predicted intensity contrasts obtained from simulations tion of the SOT spectro-polarimeter (SP) shows good agree ment, 

have also been raised and possible effects of m agnetic fields pro- even without stray light correction dDanilovic et al.ll2008h . This 

posed as explanation of the reduced contrast dUitenbroek et alJ finally should settle any remaining controversies about the RMS 

I2007I) . Detailed comparisons of the shapes of observed and granulation contrasts obtained from 3D solar convection simula- 

simulated spectral lines, determined by correlations between tions. 
Doppler velocities and temperature variation s, however, show a 

remarkable agreement dNordlund et alj |200i. In addition, there In mis paperi we initiate a search for the origin of the « miss _ 

is now excellent agreement between granulation RMS con- ing „ granulation contrast in images obserV ed with the Swedish 

trasts obtained wtth independently developed 3D MHD codes { _ meter Solar Telescope (SST) and restored with methods based 

dWedemeyer-Bohm & Rouppe van der Voort| | 2009|) . on Multi-Frame Blind Deconvolution (MFBD; iLofdahlPJOOlh. 



Recent data obtained from Hinode show higher granu- such as Joint Phase Diverse Speckle (JPDS: iPaxman et al. 



lation contrasts than obtained from most ground based so- U992) and Multi-Object MFBD (MOMFBD; Ivan Noort et al 

lar telescopes, in spite of the relatively modest aperture J2005I) . We have previously observed the effects of truncat- 

diameter of the solar optical telescope (SOT) on Hinode ing the wavefront expansion on contrasts and power spectra 

dWedemever-Bohm & Rouppe van der Voortl I2009I) . The com- when comparing data restored with MFB D based methods to 

bined PSF of the Broadband Filter Imager (BFI) and SOT, based data restored with Speckle interferometry dPaxman et alJll996t 
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iRouppe van der Voort et al.l2004l) . and this observation is the ba- 
sis of the proposed investigation. 

An important motivation for this study is Stokes measure- 
ments: Magnetic structures inside and outside sunspots have 
much more fine-structure than granulation. High spatial resolu- 
tion and low stray light is required for accurate field strength 
measurements of such structure s. Present inversion techniques, 
such as Helix (lLagg et al.ll2004t) . allow pixel-to-pixel compensa- 
tion for stray light through a parameter determined by the fits to 
the observed Stokes profiles. It is however obvious that a stray 
light PSF cannot vary significantly from one pixel to another and 
that the preferred approach is compensation for stray light in pre- 
processing and to carry out the inversions without allowance for 
stray light. 

We note that polarization signals from isolated small-scale 
structures recorded with high spatial resolution and high signal- 
to-noise (S/N) show a halo of polarized light, probably origi- 
nating from uncorrected aberrations and/or stray light. Such iso- 
lated magnetic structures provide useful "point sources" that can 
be used to validate measured PSFs in addition to measurements 
of granulation contrast at different wavelengths. 

The paper is organized as follows: We describe the 
MFBD/JPDS imaging model and reconstruction in Sect. [2] and 
propose a method for compensation of the PSF for high-order 
modes in Sect. [3] We describe simulations and tests made to val- 
idate the proposed method in Sect. [4] and apply the method to 
SST data in Sect. [5] Finally, we extend the method to include 
conventional stray light in Sect. [6] and summarize the results in 
Sect. [7] 



2. MFBD imaging model and image reconstruction 

To allow compensation for aberrations not accounted for in 
processing with MFBD based methods, we need to briefly 
review the imaging model and image reconstruction process 
used. Although the MFBD based methods include sophisticated 
schemes involving objects observed at different wavelengths 
(MOMFBD) and with known aberration differences (JPDS), the 
basic imaging model is simple and the reconstruction of the im- 
age well defined. 

In the MOMFBD implementation of Ivan Noort et alJ (120051) 
the imaging process is modeled as a space-invariant (valid for 
sufficiently small sub-fields) convolution between an unknown 
object / and a point spread function tk, the Fourier transforms of 
which are F and Tk, where the index k corresponds to a particular 
exposure. An additive Gaussian noise term »* is assumecfl The 
Fourier transform Dk of the observed image dk is then related to 
F, Tk and Nk via 

Dk = FT t +Nk. (1) 

The assumption of additive Gaussian noise leads to a maximum- 
likelihood estimate of the object and transfer functions, equiv- 
alent to the solution of a conventional non-linear least-squares 
fit problem, and corresponds to the minimization of the scalar 
quantity L, 

K 

L = J]Yj\D k -Ffk\ 2 , (2) 

H,v /;=! 



1 A maxim um-likelihood estimate with a Poisson noise model 
dPaxman et al.l|l99 2) would be more appropriate when photon noise is 
the dominating source of modelling error. The optimization process for 
the Poisson noise model is computationally much more demanding than 
for the additive Gaussian noise model. 



where u, v represents pixels in Fourier space. This equation di- 
rectly leads to an optimum estimate F of the true object F, ex- 
pressed in terms of the K estimated transfer functions 7\ and 
associated observed images dk 



^Z^/Z 1 ^ 1 



/.= ! 



(3) 



/.= ! 



dPaxman et al.ll 19921) . where IT^I 2 = tkt* k and ff is the complex 
conjugate of 7\. Note that Eq. © can be written as 



F = J] w k D k /f k 



(4) 



fc=i 



where Dk/fk are estimates of the object F based on single 
frames and w k are weights, given by 



w k = \f k \ 2 



Z^ 



(5) 



t=i 



and 2 Wk - 1 . This emphasizes that multi-frame deconvolution 
assigns weights given by the absolute squares of the estimated 
transfer functions associated with each image and at each spatial 
frequency. This is the optimum way of combining measurements 
for which S/N varies because the transfer functions vary from 
one frame to another while the RMS noise is the same for all 
measurements. In particular, zero weight is given to spatial fre- 
quencies for which the estimated transfer function pass through 
zero for an individual frame, thus avoiding division by zero as 
may happen when Dk is directly divided by Tk to obtain a re- 
stored object. Eq. (O also means that a good image Dk, corre- 
sponding to an overall large \Tk\, is given (much) higher weight 
than a poor image in the reconstruction of the object F. This 
weighting introduces a bias in the image reconstruction that is 
of importance when compensating the images for uncorrected 
high-order aberrations. 

3. Compensation for uncorrected high-order 
aberrations 

MFBD and related image reconstruction techniques do not rely 
on statistical properties of seeing, thereby allowing compensa- 
tion for individual telescope aberrations and the effects of an 
adaptive optics (AO) system in addition to seeing-induced aber- 
rations. By estimating the transfer function for each recorded 
frame, a small number of exposures can be used to restore the 
object. However, the number of aberration coefficients that can 
be determined from focused and defocused images is necessar- 
ily limited. In less than perfect seeing, this leaves a partially 
compensated PSF with enhanced wings. This corresponds to the 
well-known "halo" seen in images of point-like objects recorded 
with nigh t-time telescopes a nd partial AO compensation of see- 
ing (e.g., IConanetal.1119921) . 

3.1. Seeing statistics and residual aberrations expected 

To appreciate the importance of resid ual high-order aberrations, 
we note that iDanilovic et al.l d2008l) found that aberrations as 
small as 0.044 waves, corresponding to nearly perfect optics 
with a Strehl ratio of 93%, reduces the measured RMS con- 
trast of granulation from 8.4% to 8.1% with Hinode's spectro- 
polarimeter. Thus, within the framework of the present investi- 
gation, we clearly should consider any residual aberrations that 
decrease the Strehl ratio below 90%. 
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The effect of partially corrected aberrations on the Strehl ra- 
tio is well-known from the wo rk of Fried, Noll, Wang, Markey 
and others. We refer to Roddier ( 1999) for an overview of rele- 
vant results. The residual wavefront variance for zonal correction 
can be estimated as 



a 2 = 0.34 (D/r ) 5/3 N- 5/b , 



(6) 



where <x is in radians, D is the telescope diameter, ro is Fried's 
parameter, and N is the number of independently corrected aber- 
ration parameters. Optimum compensation with a given number 
of degrees of freedom is obtained with Karhunen-Loeve (KL) 
modeo The residual wa vefront variance , estimated from a fit to 
data plotted in Fig. 3.1 of iRoddien (I 1 9991) . is approximately 



o- 2 = 0.3 (D/r ) 5/3 AT 
The Strehl ratio R can be estimated as 
R = expC-cr 2 ). 



(1.42 



(7) 



(8) 



We shall in the following estimates assume perfect zonal correc- 
tion, which in the case of AO correction would represent an unre- 
alistically high efficiency. For the SST, with a telescope diameter 
of 0.98 m, we find that R = 0.8 1 with N = 36 aberration param- 
eters when ro = 22 cm and that R > 0.9 only when ro > 33 cm. 
With good seeing corresponding to ro = 10 cm (1" seeing) and 
very good seeing to ro = 14 cm (077 seeing), we need to com- 
pensate about 400 respectively 200 modes to reach a Strehl ratio 
of about 0.9. When ro = 20 cm (075 seeing), we need to correct 
nearly 100 modes. Thus even with excellent seeing, we should be 
concerned about the effects of uncorrected aberrations in MFBD 
processing. These estimates also imply that the problem of com- 
pensation of high-order aberrations in solar imaging is primarily 
related to near-ground seeing, associated with a large isoplanatic 
patch, since ro is typically more than a factor 3 smaller for the 
ground layer than for the high-altitude layers during day-time. 

An important limitation at short wavelengths is that ro scales 
as A 6 ^ 5 . When ro = 20 cm at A — 630 nm, ro is only 11 cm at 
390 nm. In this case we need to compensate nearly 100 aberra- 
tions at 630 nm and over 300 aberrations at 390 nm to reach a 
Strehl ratio of 0.9. Even in excellent seeing, images recorded and 
processed with conventional MFBD at 390 nm will be far from 
perfectly corrected. It is therefore not surprising that measured 
RMS contrasts of granulation with the SST show much larger 
discrepancies with 3D simulations at short wavelengths than at 
long wavelengths. 

3.2. Compensation for high-order aberrations 

In MOMFBD pro cessing of SST data , particularly from the 
CRISP instrument (Scha rmer et alj2008l) . a data set correspond- 
ing to a line scan can consist of ~ 1000 exposures in several cam- 
eras. However, an estimated object from a particular wavelength 
and polarization state is typically based on deconvolution of a 
relatively small number of images (~10), degraded by residual 
low-order aberrations from partial correction with a 37-electrode 
AO system and uncorrected high-order atmospheric aberrations. 
As long as MFBD processing compensates the residual aberra- 
tions partially corrected by the AO system, we are to some extent 
justified in ignoring the corrections made by the AO system: It is 
the residual aberrations after that AO correction that define the 
S/N of individual images and the weighting implied by Eq. (@J. 



2 This is however not true for MFBD/JPDS restoration of images ob- 
tained with AO. 



Our goal is to estimate the effect on the transfer function of 
uncompensated high-order aberrations, defined in Eq. (fTTT i. To 
do this, we will assume that these modes have amplitudes given 
by the assumption of Kolmogorov turbulence. This is the same 
basic assumption as used in Speckle Interferometry, but with the 
difference that low-order modes are compensated individually 
and for each exposed frame with MFBD. We conjecture that: 

- there will be rather small variations of the PSF from one 
frame to another because residual wavefront errors, after cor- 
rection on the order of 30 modes, depend on the accumu- 
lative effect of hundreds of modes rather than a few large- 
amplitude low-order modes. 

- a relatively small number of frames is therefore needed to 
obtain stable averages. 

Solar speckle techniques rely on an average theoretical trans- 
fer function obtained from Kolmogorov turbulence statistics and 
based on a large number of recorded frames. This corresponds 
to compensation, in an average sense, for an infinite number of 
aberration coefficients. Recent developmen ts allow solar speckle 
processin g also of AO corrected images dPuschmann & Saileri 
2006; Woger 2007) but (as far as we know) present solar speckle 
techniques do not include compensation for telescope aberra- 
tions. 



3.2.1. The proposed method 

Our goal is to develop a method that combines the advantages 
of MFBD/JPDS and speckle methods. The approach taken is to 
compensate individual frames for low-order modes and to add a 
compensation for the average effect of hundreds of high-order 
aberrations from atmospheric turbulent seeing. 

This compensation can be implemented in various ways. We 
have chosen the following simple approach that has the advan- 
tage of not involving the observed images but only their corre- 
sponding transfer functions. We note that information about the 
"true" (exact) transfer functions and the true object are encoded 
in the observed images. Ignoring the noise term in Eq. ([TJ and 
combining with Eq. ©, we obtain a relation between the esti- 
mated object F, the true object F, and the corresponding exact 
and estimated transfer functions T^ and T^ 



^=^i>^/£i^ 



k=\ 



k=\ 



We note that this is a relation of the form 



F = SF, 



(9) 



(10) 



i.e., in the form of a multiplication of the true object with a trans- 
fer function, 



A" 



k=\ 



K 



K 



z^/z^' 2 =i>4 



(id 



k=\ 



k = \ 



that is a weighted average of the ratio of the true and estimated 
transfer equations. This equation does not solve the problem un- 
less the exact transfer functions are known. However, we can 
estimate S by including known statistical properties of atmo- 
spheric high-order aberrations while using MFBD or JPDS esti- 
mates of low-order aberrations for each individual exposure. We 
propose the following estimate of S from a combination of aber- 
ration parameters determined by MFBD/JPDS processing and 
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statistical properties of atmospheric seeing: 

k k K K 

§ = Yj( T k f k / 2 \fn\ 2 ) = 2< r *> f * / Z |f * 12 (12) 

k=\ n=l k=\ k=\ 

where angular brackets, (. . .), indicate an ensemble average 
over many independent realizations for each k separately and 
where 7^ contains the same low-order aberrations, obtained with 
MFBD or JPDS processing, as 7\. 

To estimate the effect of high-order aberrations in 77, we 
add random higher-order KL modes with amplitudes given 
by Kolmogorov statistics and average the transfer equation 
over many realizations to obtain stable averages. The proposed 
method requires measurements of ro at the time of record- 
ing the data. The preferable method for accurate measure- 
ments of ro is via data from an open-loop wavefront sensor 
located before the adaptive mirror, as implemented at the SST 
dScharmer & van Werkhovenll2010l). It is also poss ible, as done 
at the Dunn telescope diameter (iMarino et al.ll2004l) . to combine 
closed-loop wavefront sensor data with the control matrix and 
output voltages to estimate ro and residual low-order aberrations. 
The quality of such measurements are to some extent limited 
by time delays and inaccuracies in the control matrix, but ex- 
perience with night-time AO systems clearly indicat es that good 
PSF c ompensation is indeed possible with such data dVeran et al.l 
fl997h . 

4. Simulations 

To investigate the feasibility of the proposed method for com- 
pensation of missing high-order aberrations in MFBD and JPDS 
processing, numerical calculations and simulations were made. 

4.1. Ideal compensation with point sources 

We us ed KL functions base d directly o n the theory of iFriedl 
(Il978h . as implemented by iDail (Il995l) . These functions are 
orthogonal on a circular aperture and statistically independent 
for Kolmogorov turbulence. For such turbulence, the variances 
depend only on ro. To produce random wavefronts following 
Kolmogorov statistics, 1001 random numbers were drawn from 
a standard normal distribution, scaled with the square root of the 
theoretical variances and used as coefficients for KL functions 
4-1004 (in decreasing variance order). The resulting wavefronts 
were scaled to different values of Fried's parameter ro. To simu- 
late the effect of partial correction with an efficient AO system, 
we reduced the amplitudes of KL coefficients 4-37 with a fac- 
tor 4 (this represents an overestimate of the actual efficiency). 
Piston (coefficient 1) is ignored because it does not contribute to 
the OTFs and tip-tilt correction (coefficients 2 & 3) was assumed 
to be perfect, corresponding to images recorded with sufficiently 
short exposure times to remove changes in image position and 
blurring during the exposure. 

The residuals of the AO-corrected aberrations were used to 
represent the estimate of the transfer equations 7\. For this calcu- 
lation we assumed 10 observed images with independently ob- 
tained wavefronts. We computed the "corrective" transfer func- 
tion S from Eq. (fTTb by using the actual high-order aberrations 
corresponding to each 7^ and then the approximate version, S , 
using the statistical averages (based on 100 realizations of the 
high-order tail) of the transfer functions, as defined in Eq. (IT2l . 
We finally compared the exact and estimated transfer functions 
S and S . This corresponds to perfect MFBD correction of the 
first 36 KL modes and perfect knowledge of ro. 



Figure Q] shows the PSFs corresponding to S and S , scaled 
such that the wings of the PSFs can be seen, for values of ro in 
the range 5-25 cm and A = 630 nm. Calculations were made 
for a 98 cm telescope diameter, corresponding to the SST, with 
critical sampling at A - 630 nm, corresponding to 07066/pixel. 
The field of view (FOV) shown is 4'.'2 x 4'.'2. As expected, the 
true PSFs show a speckled structure in the wings, but smoothed 
by the averaging effect obtained by combining ten images. The 
corresponding approximate PSFs show much less structure in 
the wings but are otherwise similar to the actual PSFs. It can 
also be seen that the effective diameter of the PSF, defined as 
that containing 90% of the energy of the PSF, increases with 
decreasing ro. When ro equals 25 cm, that diameter is approxi- 
mately 1'.' 1 , when ro is 5 cm, it increases to approximately 4'fj. 
Figure|2]shows the encircled energy as function of radius for the 
approximate and exact PSFs. It is clear that the encircled energy 
of the approximate PSF follows that of the actual PSF nearly 
exactly. In Fig. [3] we show calculated Strehl ratios of the PSFs 
as function of ro. The exact and approximate PSFs give nearly 
exactly the same Strehl ratios and also agree well with what is 
expected for perfect KL correction of 36 modes from Eqs. (0 
and (0. Finally, Fig. [4] shows radially averaged power spectra 
for the transfer functions corresponding to the exact and approx- 
imate PSFs, again showing an excellent agreement between the 
tw o at all spatial frequencies. W e also refer the reader to Fig. 12 
by iRouppe van der Voort et al . (2004), where similar effects of 
seeing on observed power spectra of penumbral fine structure 
are discussed. 



4.2. Ideal compensation with granulation images 

To investigate the effects of uncompensated high-order aberra- 
tions on granulation images, we used synth etic images calcu- 
lated from a field-free 3D MHD simulation (IStein & Nordlundl 
1 1998b ; these simulation data were kindly provided by Mats 
Carlsson. The synthetic images were calculated at a wavelength 
of 630 nm and were degraded to a resolution corresponding to 
90% of the diffraction limit of the SST while keeping the image 
scale of the original synthetic images. The images were then de- 
graded by Fourier multiplication with S based on 10 wavefronts. 
This corresponds to an MFBD estimate F based on 10 observed 
images. F was then corrected using an approximate S based on 
the same low-order aberrations as S and a high-order compen- 
sation with 1000 KL-modes according to Eq. (1121 1. The original, 
degraded and restored images are shown in Fig. [5] together with 
their RMS contrasts. Note that the degraded images differ from 
the original image in contrast but not in fine structure shown. The 
main effect of the higher-order aberrations is to add stray light, 
decreasing the RMS contrast, whereas all small-scale features of 
the original image are retained and restored to full contrast with 
the approximate PSFs. As shown in Fig. [5] the loss in contrast 
is from 14.5% to 11.8% when ro is 15 cm, or a reduction of 
the RMS contrast by nearly 20%. This corresponds to very good 
seeing conditions. In more typical seeing conditions, when ro is 
10 cm or smaller, the effects are larger. 

We conclude that whereas details of the exact and approx- 
imate PSFs certainly differ, the overall effect of the high-order 
aberrations is to add spatial stray light to the images. The pro- 



3 As discussed in Sect. 5, the increasingly larger diameter of the PSF 
in degrading seeing restricts the use of very small subfields needed with 
MFBD/JPDS methods to deal with anisoplanatism. A similar restric- 
tion applies when the images are intentionally degraded by including a 
defocused imaging channel, used with JPDS processing. 
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r = 25 cm 



20 cm 



15 cm 



10 cm 



7 cm 



5 cm 



Fig. 1: PSFs displayed in log scale. The circles mark 90% encircled energy, see also Fig. [2] Far left: Diffraction limited. Top: PSFs 
corresponding to S , i.e., true residual high-order aberrations for different ro as indicated. Bottom: Approximate PSFs, corresponding 
to S , i.e., to the method proposed. The FOV shown is 4'.'2 x 4'.'2 (64 x 64 pixels). The synthetic PSFs shown were calculated from 
10 individual noise-free frames and with perfect compensation for the 36 lowest KL-modes. 




10 15 20 

r n / 1 cm 



Fig. 2: Encircled PSF energy for different 
ro as indicated in the figure. Red: PSFs 
based on S; Blue: PSFs based on S; 
Black: diffraction limited PSF; Black 
dashed: 90% level. 



Fig. 3: Strehl ratios as a function of ro. 
Solid line: Eqs. and © with N = 37; 
Red plus (+) symbols: PSFs based on 5; 
Blue cross (x) symbols: PSFs based on S . 



Fig. 4: Power spectra (angular averages) 
of S (red) and S (blue) for ro as indicated. 



13.29% 



12.78% 



11.83% 



9.77% 



7.23% 



4.63% 




14.51% 



14.51% 



14.50% 



14.48% 



14.41% 



14.29% 



Fig. 5: Synthetic images calculated at a wavelength of 630 nm. Far left: Original image. Top: Low-pass images degraded by high- 
order aberrations, i.e., by S , corresponding to ro = 25, 20, 15, 10, 7 and 5 cm, resp. (same layout as Fig.[T). Bottom: Degraded 
images compensated by use of the method described, i.e., by S . All images are scaled between min and max of the original image 
and low-pass filtered to 90% of the SST diffraction limit. The numbers above and below the image tiles are the RMS contrasts in 
percent of the mean intensity (100 x RMS/mean). 



posed method for statistical compensation of high-order aberra- 
tions in principle should work very well. 
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4.3. Inversion tests with granulation images 

With perfect determination of the low-order modes, correction 
of the high order modes in a statistical sense gives excellent re- 
sults. However, we also need to investigate the effect of realis- 
tic errors from the low-order modes estimated with MFBD and 
JPDS techniques. 

For this test, we calculated 1000 random wavefronts cor- 
responding to Kolmogorov statistics, with reduced amplitudes 
for the first 36 KL coefficients as described above. These wave- 
fronts, scaled to different ro, were used to construct sets of de- 
graded granulation images with and without phase diversity, cor- 
responding to an added focus shift, of 1 wave peak-to-valley. 
We used synthetic images calculated from the 3D MHD code as 
true objects, compressed by a factor 2 (by cropping the Fourier 
transform of the image). The image scale was set to match that 
of SST/CRISP during the 2009 observing season, O'.'059/pixel, 
corresponding to about 12% oversampling at 630 nm. 

The degraded synthetic images were processed with the 
MOMFBD program in various ways (MFBD or JPDS, differ- 
ent numbers of realizations, different subfield sizes (256 pixels 
= 15", 128 pixels = 7/6, 80 pixels = 4'.'7), different number of 
estimated wavefront parameters, with and without added noise). 

4.3.1. Results without § compensation 

Figure [6] shows RMS intensity errors of restored images using 
different techniques (MFBD, JPDS), number of aberration pa- 
rameters (M = 35, 50 and 100) with and without S compensa- 
tion for different values of Fried's parameter ro. We emphasize 
that for all calculations made, 1000 frames were used to obtain 
the wavefronts but in panels d-f andj-1 subsets of only 10 frames 
were used to restore the images. This corresponds to processing 
of SST/CRISP data, where the broad-band channel of a MFBD 
or JPDS data set corresponds to on the order of 500-1000 im- 
ages. The simultaneously exposed narrow-band CRISP images 
are divided into typically 10-12 wavelength bins, each with 4 
polarization bins. These narrow-band images are restored indi- 
vidually, but using the aberrations (mainly) determined by the 
1000 broadband images. The upper row of plots show the results 
obtained without noise and the bottom row show the correspond- 
ing results with 0.5% Gaussian noise added to the simulated im- 
ages. 

Discussing first the results without S compensation (blue 
symbols), our reference for comparison is the results obtained 
using 1000 frames to restore the image and the first M true (ex- 
act) aberration parameters (panel a). This shows the expected 
behavior: The quality of the restored images improves with in- 
creasing value of ro and also with increasing number of esti- 
mated aberration parameters used to restore the images. When 
ro is 20 cm, the RMS error with 50 perfectly known aberration 
parameters is only 1.5%. However, in more typical seeing con- 
ditions (ro =10 cm or smaller), the corresponding RMS error 
is 2.5-8%, depending on the number of aberration parameters 
compensated and the seeing quality. 

Figure |6j5 (MFBD processing with 1000 frames) demon- 
strates that the RMS intensity error obtained using only focused 
images to estimate the aberration parameters shows quite small 
variation of the intensity error with the number of aberration 
parameters M, when M is in the range 35-100. The RMS in- 
tensity error for M — 35 MFBD calculations actually corre- 
sponds to what was obtained with the true aberration parame- 
ters for M — 100 (Panel a)! This implies that MFBD processing 
with a given number of aberration parameters leads to estimates 



of the transfer function that compensate for the effects of miss- 
ing higher-order aberration parameters (truncation of the wave- 
front). Increasing the number of estimated aberration parameters 
from 35 to 50 or 100, reduces the efficiency of this compensation 
such that the reduction of the intensity error is relatively modest. 
Quite clearly, the estimated wavefront with MFBD processing 
must be inaccurate, but such that the estimated transfer func- 
tion accurately represents the true transfer function. This conclu- 
sion is further supported by Fig. [7] which shows the variation of 
the (estimated) wavefront RMS and wavefront error for different 
number of estimated wavefront parameters when ro = 10 cm. 
Panel a (red symbols) shows that the RMS of the true wave- 
front increases by 50% when M is increased from 35 to 100. 
In contrast, the RMS of the estimated wavefront for the M first 
aberration parameters is strongly overestimated for M — 35 and 
M = 50 with MFBD processing such that the estimated wave- 
front RMS is nearly independent of M. Thus, using a low-order 
or high-order representation for the wavefront leads to nearly the 
same estimated wavefront RMS. A striking result is that MFBD 
processing of the first 35-100 corrected modes leads to residual 
wavefront errors (blue symbols) that are larger than if the wave- 
front is estimated to be exactly zero! 

When using also a defocused channel (JPDS processing), the 
wavefronts are more constrained to represent reality, leading to 
smaller wavefront errors (Fig.|7p). But this constraint also limits 
the freedom in compensating the transfer function for high-order 
aberrations, such that larger intensity errors are obtained in the 
restored images with JPDS (Fig. [(J) than with MFBD processing 

(Rg. m- 

These conclusions about differences between MFBD and 
JPDS processing are supported by the Strehl ratios calculated 
from wavefronts determined with the granulation images and 
shown in Fig. [8] With true aberration parameters (Fig. [8^) and 
without S compensation, the Strehl intensity increases grad- 
ually with the number of compensated aberration parameters 
M. However, with MFBD processing, the Strehl ratio is con- 
stant when M is in the range 35-100. The Strehl ratio achieved, 
R = 0.8, is close to what is expected from Eqs. (O and ([8]l for 
N - M - 100. The variation of the Strehl ratio with M obtained 
with JPDS processing is intermediate to that of MFBD and using 
the true aberration parameters. 

Adding 0.5% Gaussian noise to the MFBD data shows that 
the results are sensitive to noise (Figs. [6}i-i) and indicates that 
in practice the compensation effects discussed above will not be 
large and that the quality of MFBD and JPDS image restorations 
should be rather similar. The RMS intensity error improvement 
is limited to that corresponding to about 50 "true" aberration pa- 
rameters, or a Strehl ratio of 0.7. Using more than 50 aberration 
parameters to represent the wavefront does not significantly im- 
prove the quality of the restored images. 

4.3.2. Results with S compensation 

The results obtained with S compensation are shown with red 
symbols in Fig. [6] With such compensation, restorations with 
"true " low-order wavefronts allow nearly perfectly restored im- 
ages, even with 10 frames, in bad seeing (ro — 7 cm) and us- 
ing only 35 aberration parameters. This re-enforces our earlier 
conclusion, that details of the high-order part of the wavefront 
are not important as long as it has the "right" statistical proper- 
ties. In particular, the "speckles" seen in the wings of the PSF 
of Fig. Q] (upper row), have a small effect on the restored im- 
ages. The excellent result obtained with S compensation based 
on "true" wavefronts also is consistent with the results discussed 
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previously: that MFBD processing with 35 aberration parame- 
ters can produce transfer functions equivalent to those obtained 
with compensation by 100 aberration parameters, even though 
the actual wavefronts derived are demonstratedly wrong. 

With MFBD processing, the improvement from S compen- 
sation varies strongly with the number of aberration parameters, 
M. The remarkably good results obtained without S compen- 
sation for M - 35 and M = 50 constitutes a problem with S 
compensation! This is because our calculations of S do not take 
into account the compensation for high-order aberrations already 
provided by the MFBD (and JPDS) processing. This leads to 



an overcompensation of the effects of high-order aberrations for 
M - 35 and M - 50. This compensation is much smaller with 
M - 100. With JPDS processing (panels c, f, i and 1), the results 
show improved consistency. Here, S compensation with M - 50 
or M - 100 leads to RMS intensity errors that are reduced by 
a factor 2-3. The obtained RMS errors are consistently smaller 
than for MFBD processing without S compensation, when noise 
is present. 

We have here demonstrated that MFBD image reconstruc- 
tion by itself allows compensation effects by fitting transfer func- 
tions rather than wavefronts. In this sense, MFBD or JPDS pro- 
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Fig. 9: Sample wavefronts from the simulation experiments, (a)-(d): Different realizations of the Kolmogorov statistics. Within 
each sub-figure, all wavefronts are displayed with identical scaling. Left: True wavefronts; Center: JPDS estimated wavefronts; 
Right: MFBD estimated wavefronts. 1st row: 35 KL modes; 2nd row: 50 modes; 3rd row: 100 modes; 4th row: All 1001 modes. 



cessing provides more optimum image reconstruction than with 
a Shack-Hartmann based wavefront sensor. The main problem 
in the present context is that MFBD and (to a smaller extent) 
JPDS processing already introduces part of the compensation in- 
tended to be performed with S . This leads to overcompensation 
if no counter measures are taken. By using JPDS processing, the 
compensation effects are constrained and the risk of overcom- 
pensation with § is small. 

The relatively poor estimates of the wavefronts obtained with 
MFBD and JPDS processing may appear to contradict earlier 
simulations with phase div ersity methods (lLofdahl & S charmer] 
Il994t iPaxman et alJll99o) . but are a direct consequence of the 
assumed AO correction, reducing the RMS amplitudes of the 
first 35 KL aberrations with a factor of 4. This leaves residual 
wavefronts that have small or negligible amplitudes for the first 
35 modes and a correspondingly large contribution from higher- 
order modes. This leads to significant cross-talk from high-order 
aberrations, degrading the low-order wavefront estimates but ac- 
tually leading to better estimates of the transfer functions than 
expected on the basis of the number of aberration parameters in- 
cluded. Figure [9] shows the wavefronts corresponding to the first 
4 frames used to reconstruct the image with 10 frames. In the 
left column are shown (top to bottom) the true wavefronts repre- 



sented with 35, 50, 100 and 1001 KL modes. It is evident from 
the lower-left wavefront in each panel that the true wavefront has 
been stripped of its low-frequency wavefront components. It is 
also evident that the JPDS (mid column) and MFBD (right col- 
umn) wavefront estimates show large differences and that JPDS 
provides the best estimate of the true wavefront. 



4.4. Inversion tests with point sources 



We have also made inversion tests with point sources. While not 
directly relevant for restoration of solar images, these tests do 
shed further light on the MFBD and JPDS compensation effects 
discussed above. The most striking result is that both MFBD 
and JPDS processing of synthetic noise-free point source images 
lead to nearly perfectly reproduced low-order wavefronts. To un- 
derstand why granulation images and images of point sources 
lead to different wavefront estimates, we express the equation 
used to estimate the transfer function as follows: Inserting the 
expression for F in Eq. © into Eq. (0, and replacing Dk with 
FT/c, we obtain an expression for the error metric L that cor- 
responds to the scalar quantity minimized in MFBD and JPDS 
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processing: 



«,v £=1 n=l n=l 



(13) 



Ignoring the complicated expression involving the exact and es- 
timated transfer functions 7* and 7\, we emphasize that this ex- 
pression is multiplied by the power spectrum of the object, \F\ 2 . 
For a point source, this power spectrum is unity at all spatial fre- 
quencies but for solar fine structure, the power spectrum falls off 
at high spatial frequencies. Depending on the object re-imaged, 
smaller or larger weight will be given to the high spatial frequen- 
cies, obviously leading to different results as regards the derived 
wavefronts but actually quite small differences in the derived 
transfer functions. 

4.5. Comment on the use of small subfields 

The results discussed above were all obtained with 128x 128- 
pixel subfields. Calculations with 256x256-pixel subfields gave 
results that are virtually identical to those discussed. MFBD cal- 
culations with 80x80-pixel subfields gave poorer, but still ac- 
ceptable, results. However, JPDS restorations with 80x80-pixel 
subfields essentially failed. This is attributed to the relatively 
large diameter of the PSF corresponding to the defocused im- 
ages, causing problems from lack of information about the object 
outside the subfield and wrap-around effects when using FFTs 
to perform convolutions with the defocused PSF. It is thus more 
difficult to use small subfields with JPDS than with MFBD with 
the present methods. 

5. Tests with real data 

5.1. Observations 

To test the proposed S compensation method we used gran- 
ulation images recorded through a 0.44 nm FWHM 2-cavity 
filter centered at 630.26 nm and a 0.34 nm filter centered at 
538.20 nm, both used as wide-ba nd channels for th e SST/CRISP 
imaging spectro-polarimeter ( Sch armer et al. 2008). The images 
were exposed by means of a rotating chopper, set to give expo- 
sure times of 16 ms and a dark read-out time of nearly 12 ms, 
corresponding to an overall frame rate of 36 Hz. The 630 nm 
and 538 nm data discussed here were recorded between 10:45 
and 10:50 UT on 26 June 2009, during a period with reasonably 
good but strongly variable seeing. The science target for these 
observations was a small pore in AR 1023, located at approx- 
imately S22, W20, corresponding to a heliocentric distance of 
about 9 - 30° (/u = cos 9 - 0.87). This pore was also used as 
lock point for the AO system. However, in the present paper we 
discuss only field-free granulation outside this active region. The 
noise level in the recorded images was estimated from power 
spectra at high spatial frequencies and found to be about 0.9%. 

CRISP was used to repeatedly scan each of the 630.26 nm 
and 538.20 nm lines three times, using different numbers of 
wavelength positions for the two lines. These scans required a 
total of 17 s for the 630 nm line and 12 s for the 538 nm line, 
setting the time between pre-filter changes. Simultaneously, ap- 
proximately 600 wide-band images were collected at 630 nm 
and ~433 images at 538 nm during each scan. We refer to these 
17 s and 12 s sets of wide-band images as "full scan" data sets 
but will primarily discuss the results of processing 2.2 s sub-sets 
of this data to match the time-scale of seeing changes. 

The images were corrected for gain and bias and MFBD pro- 
cessed with the MOMFBD code, see Sect. IB31below. 
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Fig. 10: Variation of Fried's parameter ro with time, measured 
with the SST wide-field wavefront sensor bypassing the AO sys- 
tem (top panel) and the corresponding variation of the measured 
granulation RMS contrast of the science images (mid panel) ob- 
tained at A - 538 nm (blue curves) and at A — 630 nm (red 
curves) with the AO system in closed loop. The bottom panel 
shows the correlation between ro and the RMS contrast at the 
two wavelengths but with the contrast at A - 630 nm multiplied 
by 1.22 to bring the data onto a common curve. The dotted curve 
corresponds to the contrast of simulated data shown in Fig. 5, di- 
vided by a factor 1.85. 



10 



Scharmer, Lofdahl & al.: High-order aberration compensation with MFBD and PD image restoration techniques 



5.2. Seeing measurements 

Simultaneous measurements of Fried's parameter ro were ob- 
tained with the SST wide-field wavefront sensor (WFWFS; 
Sch armer & van Werkhoven 2010), developed as part of an ef- 
fort to characterize day-time high-altitude seeing at La Palma. 
The WFWFS is mounted on a beam that splits off light imme- 
diately before the tip-tilt and deformable mirrors, such that ro 
can be measured through the telescope and from a FOV that 
is adjacent to the science FOV, but without impact from the 
AO system. Between UT 10:00 and 10:48, processed WFWFS 
data averaged over 100 s gave estimates of ro in the range 11- 
12 cm at 500 nm, corresponding to 15-16 cm at 630 nm and 
12-13 cm at 538 nm. However, the science images, recorded 
with the AO system in closed loop, showed strong variations 
in image quality on time scales on the order of a few seconds. 

—5/3 

The WFWFS variances measured are proportional to r () (Eqs. 
(7)-(8); IScharmer & van Werk hoven 2010), such that averag- 
ing wavefront sensor data in variable seeing gives the largest 
weights to relatively poor seeing (small values of ro). However, 
MFBD image restoration gives the highest weights to good im- 
ages (large values of ro). We therefore compare data recorded 
over such short time intervals that ro can be considered con- 
stant. We re-processed the WFWFS data in blocks of 2.2 sec, 
or about 10 WFWFS CCD frames. Figure [10]shows in the up- 
per panel the variation of ro with time from UT 10:45-10:50, ob- 
tained from overlapping 2.2-sec blocks of WFWFS data. The red 
and blue curves in this panel indicate the wavelengths at which 
the science images were recorded, at 630 nm and 538 nm re- 
spectively. In the mid panel is shown the variation of the RMS 
contrast measured (see Sect. 15.4. U for the 630 nm and 538 nm 
images after flat-fielding but without MOMFBD image restora- 
tion. Comparing the two panels, there is a clear correlation be- 
tween ro and the measured RMS contrast; even small and rapid 
variations are reproduced in detail. The bottom panel shows the 
correlation between all measured values of ro and the RMS con- 
trast. In this figure, we have scaled the ro values measured with 
the WFWFS at 500 nm to 538 nm and 630 nm by assuming that 
ro is proportional to A 6 ^ 5 . We have also multiplied all RMS con- 
trasts at 630 nm by an ad hoc factor 1 .22, roughly compensating 
for the wavelength dependence of the granulation contrast. The 
plotted values cover a range of ro from 4 cm to nearly 30 cm and 
RMS contrast from 2.2% to 7.8% and show an excellent cor- 
relation. The dotted curve corresponds to variation of the RMS 
contrast with ro obtained from the simulations shown in Fig. [5] 
but divided by a factor 1.85 to fit the data. This large factor in 
part comes from the assumed perfect compensation of the first 
35 KL modes, corresponding to an AO system with 100% effi- 
ciency, and in part must be due to stray light. 

Our conclusion from the excellent correlation between mea- 
sured variations of ro with the WFWFS (bypassing the AO sys- 
tem) and the science images (recorded with the AO system in 
closed loop) is that the WFWFS indeed provides an accurate 
measure of seeing, although we cannot from this data rule out 
systematic errors, such as a scale factor error or a systematic 
bias, in the measured ro values. 



5.3. MFBD processing 

The MOMFBD software package was used to divide the im- 
ages into 128xl28-pixel overlapping subfields, which were sep- 
arately MFBD restored using the M - 100 most significant KL 
modes. The restored subfields were then mosaicked to produce 
restored images over the entire observed FOV. To investigate the 




(a) 538 nm 



(b) 630 nm 



Fig. 11: Examples of MFBD restored images based on data col- 
lected at// = 0.87 UT 10:50). The 200x400-pixel subfield shown 
was used for RMS contrast measurements. Tick marks are 1" 
apart. 



effects of strongly variable seeing on the MFBD processing, we 
processed the data in blocks corresponding to about 2.2 sec of 
data as well as in blocks corresponding to full scans (12 sec resp. 
17 sec data). To allow a comparison between the two methods of 
processing, subsets of the full scan MFBD wavefronts and ob- 
served images were used to restore the same raw images as used 
with the 2.2 sec data. 

The observed images are over-sampled by 12% at 630 nm 
but 4% under-sampled at 538 nm. In the MOMFBD code as well 
as in the S correction, under-sampling is implemented such that 
the phase over the entire pupil can be represented but the transfer 
function Tk is cropped at the Nyquist frequency. 

5.4. Results 

5.4.1. Granulation contrast measurements 

Granulation contrast was measured over the 200x400-pixel sub- 
field shown in Fig.QT] This region appears reasonably free from 
strong magnetic fields, as judged by the near absence of bright 
points and other sub-granular structure. 

The measured contrasts (blue: at 538 nm; red: at 630 nm, 
multiplied by 1.22) are shown in Fig. [12] The upper panel shows 
the contrasts measured from blocks of images recorded during 
2.2 sec intervals. The contrasts plotted in the bottom panel cor- 
respond to MFBD sets based on full scan data. The plus symbols 
refer to contrasts obtained after conventional MFBD processing 
with 100 KL aberrations, the squares to contrasts obtained after 
compensation also of high-order aberrations. We note that the 
results are quite similar for the two types of MFBD processing 
when ro is large, as expected. Data processed with 2.2 sec sets 
show a slow and gradual decrease of the RMS contrast with de- 
creasing ro after high-order compensation (squares). Inspection 
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Fig. 12: The variation of granulation RMS contrast with ro for 
MFBD restored images. Plus-symbols show the contrast af- 
ter conventional MFBD processing, squares the result after S- 
compensation. Red and blue symbols correspond to 630 nm and 
538 nm resp., with the contrasts at 630 nm multiplied by a fac- 
tor 1.22. The upper panel shows the results from processing of 
images in blocks of 2.2 sec, the bottom panel from MFBD pro- 
cessing of the entire scan to determine the wavefronts, but with 
restoration of the images and high-order compensation carried 
out in 2.2 sec. blocks. The dotted curve is identical to that shown 
in the bottom panel of Fig.fTUl 



of the restored images shows that this decrease in contrast is as- 
sociated with a small but noticeable decrease of image quality. 
For values of ro smaller than 9-10 cm, this degradation is ob- 
vious as reduced spatial resolution. A possible explanation for 
the reduced RMS contrast in poor seeing is the finite integration 
time used: A seeing layer moving at 10 ms" 1 is displaced by 
16 cm peak-to-valley (5 cm RMS) during a 16 ms exposure and 
will cause increased smearing of the wavefront with decreasing 
ro- 

For the full scan data sets obtained in strongly varying see- 
ing conditions, the trend is quite the opposite to that seen with 
the 2.2 sec data sets. For these data, the RMS contrast increases 
with decreasing ro after S compensation. This is again explained 
by the compensation effects discussed in Sect. 14.3.11 MFBD 
processing assumes that a unique object F (see Eq. [3]) is "re- 
sponsible" for all observed images of a data set, and that this 
object can be estimated with a fixed number of aberration pa- 



rameters for each of the observed images. In stable seeing, this 
works well. However, when the seeing is strongly variable, this 
causes inconsistencies. During moments of bad seeing, the miss- 
ing high-order aberrations lead to such poorly represented wave- 
fronts that the observed images are inconsistent with the im- 
ages recorded in good seeing. MFBD compensates this inconsis- 
tency by increasing the amplitudes of the low-order aberrations 
in bad seeing. When compensating for the high-order aberra- 
tions in the final reconstruction of the images, the overestimated 
wavefront RMS leads to contrast values that are too high when 
ro is small. Comparison between the restored images from the 
2.2 sec and full scan MFBD shows no differences in image qual- 
ity apart from the differences in contrast; quite clearly either way 
of MFBD processing works very well as long as the number of 
images in the data set is not too small. 



5.4.2. Comparison with 3D MHD simulations 

Based on the results of the 2.2 sec data sets, we are finally in 
a position to compare the measured granulation contrasts with 
those of simulations. The highest measured RMS contrast is 
9.8% at 538 nm and 8.5% at 630 nm (note that the contrast val- 
ues at 630 nm in Fig.[T2]are multiplied by 1.22). 

To compare measured RMS contrasts with those of the 3D 
simulations, we calculated a synthetic spectrum covering the Fe i 
630. 1 nm and 630.2 nm lines and nearby continuum, using the 
same simulation snapshot shown in Fig. [3] but for a heliocen- 
tric distance corresponding to fi - 0.87. Similarly, we included 
the weak Ci line and two stronger Fei lines for calculating a 
synthetic spectrum corresponding to the 538 nm observations. 
We multiplied the synthetic spectra obtained at each pixel with 
the transmission profiles of the CRISP pre-filters used to record 
the observed data and calculated the contrast of the granulation 
pattern. We obtained an RMS contrast of 13.4%, compared to 
13.9% at a clean nearby continuum wavelength, for the 630 nm 
data and 17.3%, compared to a continuum contrast of 17.8% 
for the 538 nm data. In comparison, the corresponding values 
at// = 1 are 14.2% and 14.6% at 630 nm and 18.2% and 18.6% 
at 538 nm. It should be emphasized that these RMS contrasts 
have been calculated without degrading the spatial resolution of 
the synthetic data to that of a 1 -m telescope, nor have any effects 
of noise on the data been included. 

The measured RMS contrasts are only 57% at 538 nm and 
63% at 630 nm of those expected, clearly demonstrating the ex- 
istence of sofar unidentified stray light sources and suggesting 
that this stray light increases at shorter wavelengths. The origin 
of this stray light will be investigated in forthcoming papers. 



5.4.3. Wavefronts 

Figure Q~3] shows the variation of the wavefront RMS with dis- 
tance from the center of the FOV, estimated from MFBD pro- 
cessing. The plot corresponds to two good data sets at each 
wavelength. Wavefront tip-tilt (KL coefficients 2 and 3) are not 
included in the calculation. As expected, the wavefront RMS 
is smallest at the center of the FOV (d = 0), where the AO is 
locked, and gradually increases away from the lock point. This 
demonstrates that in contrast to speckle techniques, MOMFBD 
techniques can be used with AO compensated images to retrieve 
wavefronts without any a priori knowledge of how anisopla- 
natism is introduced by the AO system and the Earth's atmo- 
sphere. 



12 



0.15 



© 



S 

- 0.10 



0.05 



0.00 



Scharmer, Lofdahl & al.: High-order aberration compensation with MFBD and PD image restoration techniques 

i ' ' ' ! L .' ' ' 0.8 1 ' n ' ^ I ' ' ' ' I ' ' ' ' i ' ' 



-♦♦ ! ! ! : ' ! 



+&H+T 



i 



7! 



10 



20 
d/1" 



30 



40 



Fig. 13: The variation of the wavefront RMS (in units of waves at 
500 nm) with distance d from FOV center in arc seconds. MFBD 
restorations with M — 100 corresponding to the bottom tile of 
Fig. [12] Blue: two sets of 538 nm data; Red: two sets of 630 nm 
data; Green: average of all four sets. 



Figure FBI shows the variation of the MFBD estimates of the 
wavefront RMS with ro at the center of the FOV, corresponding 
to the AO lock-point. The upper panel shows the wavefront RMS 
obtained by MFBD processing of data in 2.2 sec blocks, the bot- 
tom panel by processing full scan data. The three dotted curves 
correspond to the expected wavefront RMS with the AO sys- 
tem perfectly compensating the first 2, 9 and 16 KL aberrations 
(top to bottom, corresponding to N - 3, 10 and 17 in Eq. (|6]l). 
The top panel shows that the MFBD estimate of the wavefront 
RMS is consistent with perfect AO compensation of about 15 
KL modes when ro is approximately 10-12 cm, whereas when ro 
is larger than 20 cm, the RMS wavefront approaches the values 
expected with perfect compensation of only 10 KL modes. We 
conclude that the MFBD wavefront RMS estimates are similar, 
although somewhat higher, than what we expect from the inde- 
pendent WFWFS estimates of ro. This strongly suggests that the 
WFWFS measurements of ro are at most associated with small 
systematic errors. 

For values of ro smaller than about 9.5 cm, the MFBD esti- 
mated wavefront RMS shows a rapid transition to values equiva- 
lent to only 2 perfectly corrected KL modes. This corresponds 
to pure tip-tilt correction, or recording images with short ex- 
posures as used with the present data, but with the AO sys- 
tem not functioning at all. Such failure of solar AO systems 
to lock in poor seeing is well-known: When ro is smaller than 
the sub-aperture diameters of the Shack-Hartmann wavefront 
sensor controlling the AO system, the granulation images start 
to degrade, leading to image positions that are poorly deter- 
mined with cross-correlation techniques. The SST AO system 
uses fairly large sub-aperture diameters, 14 cm, and is more vul- 
nerable to bad seeing than other solar AO systems. The data 
in Fig. [14] suggest that the efficiency of the SST AO system 
starts to degrade when the ratio of ro to the sub-aperture diam- 
eter is about 0.7. This ratio is similar to what was found for the 
WFWFS, using sub-aperture diameters of about 9.8 cm, leading 
to strongly increased noise levels w hen ro is less than about 7 cm 
(IScharmer & van Werkhovenll2010l) . 

The bottom panel in Fig. [14] based on full-scan data recorded 
in strongly variable seeing, shows a variation of wavefront RMS 
with ro that is quite similar to that of the upper panel. However, 
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Fig. 14: RMS wavefronts vs ro. Red and blue circles: RMS wave- 
fronts as estimated with MFBD using M — 100 modes and 
2.2 sec data blocks (top tile), and 12-17 sec data blocks (bot- 
tom tile). Black dotted curves: Wavefront RMS expected from 
perfect correction of the first N —\ modes for N = 3, 10, 17 (top 
to bottom) KL modes. 



the wavefront RMS, in particular during bad seeing, is systemati- 
cally higher for the full-scan data. This is consistent with the sys- 
tematic differences in granulation contrasts found for these data, 
discussed in the previous section: In strongly varying seeing the 
limited number of aberrations used to model the wavefront be- 
comes a problem in particular during moments of poor seeing. 
The MFBD compensates this by increasing the RMS of the low- 
order aberrations. When combined with S -compensation, this 
leads to overcompensation. At the same time, MFBD optimiza- 
tion is dominated by the best frames, because of the weighting 
of the transfer functions used to restore the object, see Eq. [4] 
This reduces the contribution of the poorer and overcompensated 
frames in the restored image and thus leads to a reduction of the 
overcompensation effects when images recorded in variable see- 
ing are combined and restored. 



6. Stray-Light compensation 

An in-depth discussion of stray light measurements and com- 
pensation is beyond the scope of the present paper. However, we 
note that the primary focal plane used to calibrate the control ma- 
trix of an AO system can also be used to aid stray light calibra- 
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tions of any following optics and instrumentation. Such (partial) 
stray light calibrations can be made by locking the AO system 
on the stray light target and recording images with the science 
camera at several focus positions. Processing these images with 
JPDS techniques, allows the effects of any residual aberrations 
to be separated from those of stray light, making such measure- 
ments robust. 

Spatial stray light can be simply modeled as a convolution: 



D k =FT s T k + N k , 



(14) 



where T s is the Fourier transform of the stray light PSF, t s . The 
residual aberrations can be identified and determined by record- 
ing pinhole images at different focus locations and the corre- 
sponding transfer function 7* can be calculated and compen- 
sated for with phase diversity techniques. The stray light transfer 
function T s is assumed to be invariant with respect to the focus 
position. It is then evident that phase diversity image restoration 
will produce a restored "object" that is the product of F and T s 
and that T s can therefore only be determined if the object F is 
known. For example, using a circular pinhole as object, corre- 
sponds to a binary object that is of constant intensity within a 
radius r and zero outside that radius. Recording images of pin- 
holes of different diameters should allow the stray light PSF to 
be determined with confidence. Procedures for such stray light 
measurements and compensation will be described in a forth- 
coming paperQ 

Of relevance to the present paper is that the stray light PSF, 
as well as S in bad seeing, have wings that can extend well be- 
yond several arc sec. Proper compensation for such a broad PSF 
during MFBD/JPDS processing would cause problems with the 
small sub-fields needed to deal with anisoplanatism. The calcu- 
lation of S involves atmospheric transfer functions modified by 
high-order aberrations but does not involve the observed images 
dk- A reasonable expectation is that S should show relatively 
small and smooth variations over a large FOV. This would allow 
S T s , averaged over a large fraction over the FOV, to compensate 
for the combined effects of high-order atmospheric aberrations 
and stray light using only a few deconvolutions, independent of 
the sub-fielding used in previous MFBD or JPDS processing. 
The restored images can then be properly compensated for stray 
light originating from the broad wings of the PSF corresponding 
to ST ' s . 

7. Conclusions 

The development of techniques for restoration of solar im- 
ages based on multi-fr ame blind deconcolution (M FBD) such 
as phase diversity (PD; Lofdahl & Scharmei 1994b . joint phase 
diversity speckle (JPDS: IPaxman et aDll992|) a nd multi-object 
multi-frame blind de convolution (MOMFBD; lLofdahll l2002t 
Ivan Noort et al.l 120051) is of major importance to present and 
future broad-band imaging, imaging spectrometry and imaging 
spectropolarimetry with ground-based solar telescopes. An im- 
portant advantage of these techniques is that they do not rely on 



4 We remark that such a calibration will also include the detector 
MTF. This MTF is primarily the combination of two effects: The first 
is from the integration over the finite size of the pixels and the second 
is fro m charge diffusion in the silicon substrate of the CCD detector 
(e.g., IStevens & Lavinej[l994h . With back-illuminated thinned CCDs 
operating at near infrared wavelengths, as used with CRISP, there is a 
third stray light source from light that is first transmitted through the 
thin silicon layer of the CCD and then scattered back into the beam. 
Compensation for such stray light needs special calibration and com- 
pensation techniques (de la Cruz Rodriguez et al., in preparation). 



statistical properties of atmospheric seeing, including anisopla- 
natism introduced by the Earth's atmosphere and a conventional 
or multi-conjugate AO system. These techniqu es, implemented 
in the MOMFBD C++ code developed by Ivan Noort et al.l 
(120051) . have been used extensively to process data sets recorded 
with the SST This has resulted in stunning time sequences with 
near-diffraction limited resolution. The remarkably stable image 
quality achieved when many exposed frames are processed as a 
single data set strongly indicates that this code is working well. 
The present work represents the first systematic evaluation of its 
performance. 

Our simulations confirm and explain the excellent perfor- 
mance of the code with images recorded using a low-order (37 
electrode) AO system. In particular, we demonstrate that even 
though MFBD processing of AO compensated images leads to 
poorly determined wavefronts, this is of secondary importance. 
The transfer functions, needed to restore images, are actually 
better determined than expected from the number of KL modes 
used to represent the wavefront. A surprising result is that the 
addition of a phase diversity (defocused) channel not necessarily 
improves the result when many images are processed simultane- 
ously: adding such a channel constrains the wavefronts to more 
closely represent reality, but the obtained transfer functions may 
actually be more accurate without the defocused channel. From 
this perspective, MFBD processing rather than JPDS process- 
ing is preferrable. However, when applying § -compensation to 
the restored images, the more accurate wavefronts of JPDS con- 
stitute an advantage. A disadvantage of JPDS processing with 
an added diversity channel is that defocused images are fuzzier, 
causing "leakage" of (stray-)light across subfield boundaries and 
wrapping around when using FFTs to perform convolutions. 
This degrades the quality of the restorations or even cause the 
inversions to fail when the subfields are small. We can thus not 
simply conclude that the MFBD or JPDS method of processing 
in general is preferrable, but always recording also defocused 
images is a good strategy since this allows either mode of pro- 
cessing the data. 

A drawback of MFBD and JPDS techniques is that trunca- 
tion of the wavefront expansion is necessary. This leaves a tail of 
uncompensated high-order aberrations and under-compensation 
of the restored images. Based on restorations of images recorded 
in strongly variable seeing, we have found that MFBD process- 
ing compensates for the missing high-order aberrations by in- 
creasing the amplitudes of the low-order aberrations. This com- 
pensation effect is particularly large during moments of bad see- 
ing but the corresponding exposed frames are given relatively 
small weight in the restored images. This explains the appar- 
ent stable image quality of movies made from MFBD processed 
data. 

We have proposed a method for further compensation 
of missing high-order aberrations, based on statistics of 
Kolmogorov turbulence. Our simulations show that such com- 
pensation in principle should allow perfect image restorations 
in good to excellent seeing conditions. A major advantage of 
the method is that very good results can be obtained even when 
only a small number of frames are used to restore an object. 
However, the high-order mode compensations already intro- 
duced by MFBD or JPDS image processing limits the overall 
improvement factor to about 2-3, which nevertheless should rep- 
resent an important improvement of the photometric quality of 
the restored images. A problem is MFBD processing of data 
sets recorded during strongly variable seeing. In such conditions, 
the use of a truncated wavefront leads to overcompensation for 
the images recorded during moments of poor seeing such that 
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the contrast of the restored images increases with decreasing ro 
when applying S -compensation. By breaking up the processing 
in smaller data sets, these overcompensation effects are strongly 
reduced. 

Using images recorded with the SST in 1" seeing, we have 
calculated the RMS contrast of solar granulation. The obtained 
"raw" contrast values range from 7.5% at 538 nm and 6.5% at 
630 nm for the best images, observed away from sun center (at 
cos 9 - /j. — 0.87) and through pre-filters that contain a few mod- 
erately strong Fe i lines. Due to the previous compensation by 
the AO system, the measured RMS contrast increases by only 
a fairly small amount to 9% and 7.5% respectively for MFBD 
processed images using 100 KL aberrations. Simultaneous mea- 
surements of Fried's parameter ro were used to aid compensa- 
tion for high-order aberrations and yields an increase of the RMS 
contrast to close to 10% at 538 nm and 8.5% at 630 nm. The cor- 
responding RMS contrast of synthetic images obtained from 3D 
MHD simulations, taking into account the strongest lines within 
the filters used and an inclination angle of 6 — 30°, is 17.4% 
at 538 nm and 13.3% at 630 nm. This implies a remaining dis- 
crepancy of 37-43% of the true contrast, the major part of which 
must be from stray light in the Earth's atmosphere, the telescope, 
the subsequent re-imaging optics, the finite (16 ms) integration 
time used to record the images and/or anisoplanatism, scintil- 
lation and image scale variations from high-altitude seeing. The 
origin of this missing stray light will be the topic of forth-coming 
papers. 
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